#ifndef __MATH_H
#define __MATH_H
typedef struct
{
    int quot, rem;
} div_t;
div_t divmod (int numer, int denom)
{
    div_t result;
    result.quot = numer / denom;
    result.rem = numer % denom;
    return result;
}

#define E 2.7182818284590452353602874713526624977572470936999595749
#define EPSILON 0.0000000000000000001
#define ONE 1.0
#define PI 3.14159265358979323846264338327950288419716939937510
#define HALF_PI 1.570796326794896619231321691639751442098584699687552910487
#define ZERO 0.0
#define E_INVERTED 0.3678794411714423215955237701614608674458111310317678
#define HALF 0.5
#define LOG10_INV 0.434294481903251827651128918916605082294397005803666566114
#define MAX_ITERATIONS 100
#define QUARTER_PI 0.785398163397448309615660845819875721049292349843776455243
double abs(double x) {
    if (x <= ZERO) {
        return -x;
    }
    return x;
}
double exp(double x) {
    int count = 0;
    while (x > E) {
        x--;
        count++;
    }

    while (x < ZERO) {
        x++;
        count--;
    }

    x--;

    if (x == 0) {
        return count;
    }

    double result = ZERO;
    double iteration = 0;
    double y = ONE;
    double cacheResult = result - ONE;
    while (cacheResult != result && iteration < MAX_ITERATIONS) {
        iteration++;
        cacheResult = result;
        y *= -x;
        result += y / iteration;
    }

    return count - result;
}



double asin(double x) {
    if (x > ONE || x < -ONE) {
        return NULL;
    }

    if (x == ZERO) {
        return ZERO;
    } else if (x == ONE) {
        return HALF_PI;
    } else {
        x *= x;
        double y = ONE;
        double cachedY = y - ONE;
        for (int i = 1; cachedY != y && i < MAX_ITERATIONS; i++) {
            cachedY = y;
            y += x * (2 * i - 1) / (2 * i);
        }
        return y;
    }
}


double cosh(double x) {
    double y = exp(x);
    double yy = ONE / y;
    return (y + yy) * HALF;
}


double log(double x) {
    if (x <= ZERO) {
        return NULL;
    }

    int count = 0;
    while (x > E) {
        x /= E;
        count++;
    }

    double result = ZERO;
    double iteration = 0;
    double y = ONE;
    double cacheResult = result - ONE;
    while (cacheResult != result && iteration < MAX_ITERATIONS) {
        iteration++;
        cacheResult = result;
        y *= x;
        result += y / iteration;
    }

    return -count - result;
}

double log10(double x) {
    if (x <= ZERO) {
        return NULL;
    }

    int count = 0;
    while (x > E) {
        x /= E;
        count++;
    }

    double result = ZERO;
    double iteration = 0;
    double y = ONE;
    double cacheResult = result - ONE;
    while (cacheResult != result && iteration < MAX_ITERATIONS) {
        iteration++;
        cacheResult = result;
        y *= x;
        result += y / iteration;
    }

    return -count - result * LOG10_INV;
}
int quickpow(int x,int y){
    if(y==0) return 1;
    if(y==1) return x;
    int result = 1;
    while (y > 0) {
        if (y % 2 == 0) {
            x *= x;
            y /= 2;
        } else {
            result *= x;
            y--;
        }
    }
    return result;
}
double pow(double x, double y) {
    if (x == ZERO && y == ZERO) {
        return NULL;
    }

    if (x == ZERO) {
        return ZERO;
    }

    if (y == ZERO) {
        return ONE;
    }

    if (x == ONE) {
        return y;
    }

    if (y == ONE) {
        return x;
    }

    if (x > E && y > E) {
        return exp(y * log(x));
    }

    if (x < E && y < E) {
        return exp(-y * log(x));
    }

    double result = ZERO;
    double iteration = 0;
    double yInv = ONE / y;
    double cacheResult = result - ONE;
    while (cacheResult != result && iteration < MAX_ITERATIONS) {
        iteration++;
        cacheResult = result;
        result += yInv * pow(x, iteration);
    }

    return result;
}

double sinh(double x) {
    double y = exp(x);
    double yy = ONE / y;
    return (y - yy) * HALF;
}

double sqrt(double x) {
    if (x < ZERO) {
        return NULL;
    }

    if (x == ZERO) {
        return ZERO;
    }

    double xhalf = x;
    int i = 0;
    double epsilon = EPSILON;

    while (abs(xhalf * xhalf - x) > epsilon) {
        i++;
        xhalf = (xhalf + x) / 2;
    }

    return xhalf;
}

double tanh(double x) {
    double y = exp(x);
    double yy = ONE / y;
    return (y - yy) / (y + yy);
}

double acos(double x) {
    if (x == ZERO) {
        return HALF_PI;
    } else if (x == ONE) {
        return ZERO;
    } else if (x < ZERO) {
        return PI - acos(-x);
    } else {
        return HALF_PI - asin(x);
    }
}
double atan(double x) {
    if (x == ZERO) {
        return ZERO;
    } else if (x == ONE) {
        return QUARTER_PI;
    } else {
        return asin(x / sqrt(ONE + x * x));
    }
}
#endif